home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / SUPER.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-06-26  |  29.7 KB  |  1,648 lines  |  [TEXT/CWIE]
  1. /****************************************************************************
  2. *                   super.c
  3. *
  4. *  This module implements functions that manipulate superellipsoids.
  5. *
  6. *  Original code written by Alexander Enzmann.
  7. *  Adaption to POV-Ray by Dieter Bayer [DB].
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. /****************************************************************************
  28. *
  29. *  Explanation:
  30. *
  31. *    Superellipsoids are defined by the implicit equation
  32. *
  33. *      f(x,y,z) = (|x|^(2/e) + |y|^(2/e))^(e/n) + |z|^(2/n) - 1
  34. *
  35. *    Where e is the east/west exponent and n is the north/south exponent.
  36. *
  37. *    NOTE: The exponents are precomputed and stored in the Power element.
  38. *
  39. *    NOTE: Values of e and n that are close to degenerate (e.g.,
  40. *          below around 0.1) appear to give the root solver fits.
  41. *          Not sure quite where the problem lays just yet.
  42. *
  43. *  Syntax:
  44. *
  45. *    superellipsoid { <e, n> }
  46. *
  47. *
  48. *  ---
  49. *
  50. *  Oct 1994 : Creation.
  51. *
  52. *****************************************************************************/
  53.  
  54. #include "frame.h"
  55. #include "povray.h"
  56. #include "vector.h"
  57. #include "povproto.h"
  58. #include "bbox.h"
  59. #include "matrices.h"
  60. #include "objects.h"
  61. #include "super.h"
  62.  
  63.  
  64.  
  65. /*****************************************************************************
  66. * Local preprocessor defines
  67. ******************************************************************************/
  68.  
  69. /* Minimal intersection depth for a valid intersection. */
  70.  
  71. #define DEPTH_TOLERANCE 1.0e-4
  72.  
  73. /* If |x| < ZERO_TOLERANCE it is regarded to be 0. */
  74.  
  75. #define ZERO_TOLERANCE 1.0e-10
  76.  
  77. /* This is not the signum function because SGNX(0) is 1. */
  78.  
  79. #define SGNX(x) (((x) >= 0.0) ? 1.0 : -1.0)
  80.  
  81. #define MIN_VALUE -1.0
  82. #define MAX_VALUE  1.0
  83.  
  84. #define MAX_ITERATIONS 20
  85.  
  86. #define PLANECOUNT 9
  87.  
  88.  
  89.  
  90. /*****************************************************************************
  91. * Local typedefs
  92. ******************************************************************************/
  93.  
  94.  
  95.  
  96. /*****************************************************************************
  97. * Static functions
  98. ******************************************************************************/
  99.  
  100. static int intersect_superellipsoid PARAMS((RAY *Ray, SUPERELLIPSOID *Superellipsoid, ISTACK *Depth_Stack));
  101. static int intersect_box PARAMS((VECTOR P, VECTOR D, DBL *dmin, DBL *dmax));
  102. static DBL power PARAMS((DBL x, DBL e));
  103. static DBL evaluate_superellipsoid PARAMS((VECTOR P, SUPERELLIPSOID *Superellipsoid));
  104. static int compdists PARAMS((CONST void *in_a, CONST void *in_b));
  105.  
  106. static int find_ray_plane_points PARAMS((VECTOR P, VECTOR D, int cnt, DBL *dists, DBL mindist, DBL maxdist));
  107. static void solve_hit1 PARAMS((SUPERELLIPSOID *Super, DBL v0, VECTOR tP0, DBL v1, VECTOR tP1, VECTOR P));
  108. static int check_hit2 PARAMS((SUPERELLIPSOID *Super, VECTOR P, VECTOR D, DBL t0, VECTOR P0, DBL v0, DBL t1, DBL *t, VECTOR Q));
  109. static int insert_hit PARAMS((OBJECT *Object, RAY *Ray, DBL Depth, ISTACK *Depth_Stack));
  110.  
  111. static int   All_Superellipsoid_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  112. static int   Inside_Superellipsoid PARAMS((VECTOR point, OBJECT *Object));
  113. static void  Superellipsoid_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  114. static void  *Copy_Superellipsoid PARAMS((OBJECT *Object));
  115. static void  Translate_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  116. static void  Rotate_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  117. static void  Scale_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  118. static void  Transform_Superellipsoid PARAMS((OBJECT *Object, TRANSFORM *Trans));
  119. static void  Invert_Superellipsoid PARAMS((OBJECT *Object));
  120. static void  Destroy_Superellipsoid PARAMS((OBJECT *Object));
  121.  
  122. /*****************************************************************************
  123. * Local variables
  124. ******************************************************************************/
  125.  
  126. static METHODS Superellipsoid_Methods =
  127. {
  128.   All_Superellipsoid_Intersections,
  129.   Inside_Superellipsoid, Superellipsoid_Normal,
  130.   Copy_Superellipsoid,
  131.   Translate_Superellipsoid, Rotate_Superellipsoid,
  132.   Scale_Superellipsoid, Transform_Superellipsoid, Invert_Superellipsoid, Destroy_Superellipsoid
  133. };
  134.  
  135. static DBL planes[PLANECOUNT][4] =
  136.   {{1, 1, 0, 0}, {1,-1, 0, 0},
  137.    {1, 0, 1, 0}, {1, 0,-1, 0},
  138.    {0, 1, 1, 0}, {0, 1,-1, 0},
  139.    {1, 0, 0, 0}, 
  140.    {0, 1, 0, 0}, 
  141.    {0, 0, 1, 0}};
  142.  
  143.  
  144.  
  145. /*****************************************************************************
  146. *
  147. * FUNCTION
  148. *
  149. *   All_Superellipsoid_Intersections
  150. *
  151. * INPUT
  152. *
  153. *   Object      - Object
  154. *   Ray         - Ray
  155. *   Depth_Stack - Intersection stack
  156. *   
  157. * OUTPUT
  158. *
  159. *   Depth_Stack
  160. *   
  161. * RETURNS
  162. *
  163. *   int - TRUE, if a intersection was found
  164. *   
  165. * AUTHOR
  166. *
  167. *   Dieter Bayer
  168. *   
  169. * DESCRIPTION
  170. *
  171. *   Determine ray/superellipsoid intersection and clip intersection found.
  172. *
  173. * CHANGES
  174. *
  175. *   Oct 1994 : Creation.
  176. *
  177. ******************************************************************************/
  178.  
  179. static int All_Superellipsoid_Intersections(Object, Ray, Depth_Stack)
  180. OBJECT *Object;
  181. RAY *Ray;
  182. ISTACK *Depth_Stack;
  183. {
  184.   Increase_Counter(stats[Ray_Superellipsoid_Tests]);
  185.  
  186.   if (intersect_superellipsoid(Ray, (SUPERELLIPSOID *)Object, Depth_Stack))
  187.   {
  188.     Increase_Counter(stats[Ray_Superellipsoid_Tests_Succeeded]);
  189.  
  190.     return(TRUE);
  191.   }
  192.   else
  193.   {
  194.     return(FALSE);
  195.   }
  196. }
  197.  
  198.  
  199.  
  200. /*****************************************************************************
  201. *
  202. * FUNCTION
  203. *
  204. *   intersect_superellipsoid
  205. *
  206. * INPUT
  207. *
  208. *   Ray            - Ray
  209. *   Superellipsoid - Superellipsoid
  210. *   Depth_Stack    - Depth stack
  211. *   
  212. * OUTPUT
  213. *
  214. *   Intersection
  215. *   
  216. * RETURNS
  217. *
  218. *   int - TRUE if intersections were found.
  219. *   
  220. * AUTHOR
  221. *
  222. *   Alexander Enzmann, Dieter Bayer
  223. *   
  224. * DESCRIPTION
  225. *
  226. *   Determine ray/superellipsoid intersection.
  227. *
  228. * CHANGES
  229. *
  230. *   Oct 1994 : Creation.
  231. *
  232. ******************************************************************************/
  233.  
  234. static int intersect_superellipsoid(Ray, Superellipsoid, Depth_Stack)
  235. RAY *Ray;
  236. SUPERELLIPSOID *Superellipsoid;
  237. ISTACK *Depth_Stack;
  238. {
  239.   int i, cnt, Found = FALSE;
  240.   DBL dists[PLANECOUNT+2];
  241.   DBL t, t1, t2, v0, v1, len;
  242.   VECTOR P, D, P0, P1, P2, P3;
  243.  
  244.   /* Transform the ray into the superellipsoid space. */
  245.  
  246.   MInvTransPoint(P, Ray->Initial, Superellipsoid->Trans);
  247.  
  248.   MInvTransDirection(D, Ray->Direction, Superellipsoid->Trans);
  249.  
  250.   VLength(len, D);
  251.  
  252.   VInverseScaleEq(D, len);
  253.  
  254.   /* Intersect superellipsoid's bounding box. */
  255.  
  256.   if (!intersect_box(P, D, &t1, &t2))
  257.   {
  258.     return(FALSE);
  259.   }
  260.  
  261.   /* Test if superellipsoid lies 'behind' the ray origin. */
  262.  
  263.   if (t2 < DEPTH_TOLERANCE)
  264.   {
  265.     return(FALSE);
  266.   }
  267.  
  268.   cnt = 0;
  269.  
  270.   if (t1 < DEPTH_TOLERANCE)
  271.   {
  272.     t1 = DEPTH_TOLERANCE;
  273.   }
  274.  
  275.   dists[cnt++] = t1;
  276.   dists[cnt++] = t2;
  277.  
  278.   /* Intersect ray with planes cutting superellipsoids in pieces. */
  279.  
  280.   cnt = find_ray_plane_points(P, D, cnt, dists, t1, t2);
  281.  
  282.   if (cnt <= 1)
  283.   {
  284.     return(FALSE);
  285.   }
  286.  
  287.   VEvaluateRay(P0, P, dists[0], D)
  288.  
  289.   v0 = evaluate_superellipsoid(P0, Superellipsoid);
  290.  
  291.   if (fabs(v0) < ZERO_TOLERANCE)
  292.   {
  293.     if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[0] / len, Depth_Stack))
  294.     {
  295.       if (Superellipsoid->Type & IS_CHILD_OBJECT)
  296.       {
  297.         Found = TRUE;
  298.       }
  299.       else
  300.       {
  301.         return(TRUE);
  302.       }
  303.     }
  304.   }
  305.  
  306.   for (i = 1; i < cnt; i++)
  307.   {
  308.     VEvaluateRay(P1, P, dists[i], D)
  309.  
  310.     v1 = evaluate_superellipsoid(P1, Superellipsoid);
  311.  
  312.     if (fabs(v1) < ZERO_TOLERANCE)
  313.     {
  314.       if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[i] / len, Depth_Stack))
  315.       {
  316.         if (Superellipsoid->Type & IS_CHILD_OBJECT)
  317.         {
  318.           Found = TRUE;
  319.         }
  320.         else
  321.         {
  322.           return(TRUE);
  323.         }
  324.       }
  325.     }
  326.     else
  327.     {
  328.       if (v0 * v1 < 0.0)
  329.       {
  330.         /* Opposite signs, there must be a root between */
  331.  
  332.         solve_hit1(Superellipsoid, v0, P0, v1, P1, P2);
  333.  
  334.         VSub(P3, P2, P);
  335.  
  336.         VLength(t, P3);
  337.  
  338.         if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
  339.         {
  340.           if (Superellipsoid->Type & IS_CHILD_OBJECT)
  341.           {
  342.             Found = TRUE;
  343.           }
  344.           else
  345.           {
  346.             return(TRUE);
  347.           }
  348.         }
  349.       }
  350.       else
  351.       {
  352.         /* 
  353.          * Although there was no sign change, we may actually be approaching
  354.          * the surface. In this case, we are being fooled by the shape of the
  355.          * surface into thinking there isn't a root between sample points. 
  356.          */
  357.  
  358.         if (check_hit2(Superellipsoid, P, D, dists[i-1], P0, v0, dists[i], &t, P2))
  359.         {
  360.           if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
  361.           {
  362.             if (Superellipsoid->Type & IS_CHILD_OBJECT)
  363.             {
  364.               Found = TRUE;
  365.             }
  366.             else
  367.             {
  368.               return(TRUE);
  369.             }
  370.           }
  371.           else
  372.           {
  373.             break;
  374.           }
  375.         }
  376.       }
  377.     }
  378.  
  379.     v0 = v1;
  380.  
  381.     Assign_Vector(P0, P1);
  382.   }
  383.  
  384.   return(Found);
  385. }
  386.  
  387.  
  388.  
  389. /*****************************************************************************
  390. *
  391. * FUNCTION
  392. *
  393. *   Inside_Superellipsoid
  394. *
  395. * INPUT
  396. *
  397. *   IPoint - Intersection point
  398. *   Object - Object
  399. *   
  400. * OUTPUT
  401. *   
  402. * RETURNS
  403. *
  404. *   int - TRUE if inside
  405. *   
  406. * AUTHOR
  407. *
  408. *   Dieter Bayer
  409. *   
  410. * DESCRIPTION
  411. *
  412. *   Test if a point lies inside the superellipsoid.
  413. *
  414. * CHANGES
  415. *
  416. *   Oct 1994 : Creation.
  417. *
  418. ******************************************************************************/
  419.  
  420. static int Inside_Superellipsoid(IPoint, Object)
  421. VECTOR IPoint;
  422. OBJECT *Object;
  423. {
  424.   DBL val;
  425.   VECTOR P;
  426.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  427.  
  428.   /* Transform the point into the superellipsoid space. */
  429.  
  430.   MInvTransPoint(P, IPoint, Superellipsoid->Trans);
  431.  
  432.   val = evaluate_superellipsoid(P, Superellipsoid);
  433.  
  434.   if (val < EPSILON)
  435.   {
  436.     return(!Test_Flag(Superellipsoid, INVERTED_FLAG));
  437.   }
  438.   else
  439.   {
  440.     return(Test_Flag(Superellipsoid, INVERTED_FLAG));
  441.   }
  442. }
  443.  
  444.  
  445.  
  446. /*****************************************************************************
  447. *
  448. * FUNCTION
  449. *
  450. *   Superellipsoid_Normal
  451. *
  452. * INPUT
  453. *
  454. *   Result - Normal vector
  455. *   Object - Object
  456. *   Inter  - Intersection found
  457. *   
  458. * OUTPUT
  459. *
  460. *   Result
  461. *   
  462. * RETURNS
  463. *   
  464. * AUTHOR
  465. *
  466. *   Dieter Bayer
  467. *   
  468. * DESCRIPTION
  469. *
  470. *   Calculate the normal of the superellipsoid in a given point.
  471. *
  472. * CHANGES
  473. *
  474. *   Oct 1994 : Creation.
  475. *
  476. ******************************************************************************/
  477.  
  478. static void Superellipsoid_Normal(Result, Object, Inter)
  479. OBJECT *Object;
  480. VECTOR Result;
  481. INTERSECTION *Inter;
  482. {
  483.   DBL k, x, y, z;
  484.   VECTOR P, N, E;
  485.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  486.  
  487.   /* Transform the point into the superellipsoid space. */
  488.  
  489.   MInvTransPoint(P, Inter->IPoint, Superellipsoid->Trans);
  490.  
  491.   Assign_Vector(E, Superellipsoid->Power);
  492.  
  493.   x = fabs(P[X]);
  494.   y = fabs(P[Y]);
  495.   z = fabs(P[Z]);
  496.  
  497.   k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
  498.  
  499.   N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
  500.   N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
  501.   N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
  502.  
  503.   /* Transform the normalt out of the superellipsoid space. */
  504.  
  505.   MTransNormal(Result, N, Superellipsoid->Trans);
  506.  
  507.   VNormalize(Result, Result);
  508. }
  509.  
  510.  
  511.  
  512. /*****************************************************************************
  513. *
  514. * FUNCTION
  515. *
  516. *   Translate_Superellipsoid
  517. *
  518. * INPUT
  519. *
  520. *   Object - Object
  521. *   Vector - Translation vector
  522. *
  523. * OUTPUT
  524. *
  525. *   Object
  526. *
  527. * RETURNS
  528. *
  529. * AUTHOR
  530. *
  531. *   Dieter Bayer
  532. *
  533. * DESCRIPTION
  534. *
  535. *   Translate a superellipsoid.
  536. *
  537. * CHANGES
  538. *
  539. *   Oct 1994 : Creation.
  540. *
  541. ******************************************************************************/
  542.  
  543. static void Translate_Superellipsoid(Object, Vector, Trans)
  544. OBJECT *Object;
  545. VECTOR Vector;
  546. TRANSFORM *Trans;
  547. {
  548.   Transform_Superellipsoid(Object, Trans);
  549. }
  550.  
  551.  
  552.  
  553. /*****************************************************************************
  554. *
  555. * FUNCTION
  556. *
  557. *   Rotate_Superellipsoid
  558. *
  559. * INPUT
  560. *
  561. *   Object - Object
  562. *   Vector - Rotation vector
  563. *   
  564. * OUTPUT
  565. *
  566. *   Object
  567. *   
  568. * RETURNS
  569. *   
  570. * AUTHOR
  571. *
  572. *   Dieter Bayer
  573. *   
  574. * DESCRIPTION
  575. *
  576. *   Rotate a superellipsoid.
  577. *
  578. * CHANGES
  579. *
  580. *   Oct 1994 : Creation.
  581. *
  582. ******************************************************************************/
  583.  
  584. static void Rotate_Superellipsoid(Object, Vector, Trans)
  585. OBJECT *Object;
  586. VECTOR Vector;
  587. TRANSFORM *Trans;
  588. {
  589.   Transform_Superellipsoid(Object, Trans);
  590. }
  591.  
  592.  
  593.  
  594. /*****************************************************************************
  595. *
  596. * FUNCTION
  597. *
  598. *   Scale_Superellipsoid
  599. *
  600. * INPUT
  601. *
  602. *   Object - Object
  603. *   Vector - Scaling vector
  604. *   
  605. * OUTPUT
  606. *
  607. *   Object
  608. *   
  609. * RETURNS
  610. *   
  611. * AUTHOR
  612. *
  613. *   Dieter Bayer
  614. *   
  615. * DESCRIPTION
  616. *
  617. *   Scale a superellipsoid.
  618. *
  619. * CHANGES
  620. *
  621. *   Oct 1994 : Creation.
  622. *
  623. ******************************************************************************/
  624.  
  625. static void Scale_Superellipsoid(Object, Vector, Trans)
  626. OBJECT *Object;
  627. VECTOR Vector;
  628. TRANSFORM *Trans;
  629. {
  630.   Transform_Superellipsoid(Object, Trans);
  631. }
  632.  
  633.  
  634.  
  635. /*****************************************************************************
  636. *
  637. * FUNCTION
  638. *
  639. *   Transform_Superellipsoid
  640. *
  641. * INPUT
  642. *
  643. *   Object - Object
  644. *   Trans  - Transformation to apply
  645. *   
  646. * OUTPUT
  647. *
  648. *   Object
  649. *   
  650. * RETURNS
  651. *   
  652. * AUTHOR
  653. *
  654. *   Dieter Bayer
  655. *   
  656. * DESCRIPTION
  657. *
  658. *   Transform a superellipsoid and recalculate its bounding box.
  659. *
  660. * CHANGES
  661. *
  662. *   Oct 1994 : Creation.
  663. *
  664. ******************************************************************************/
  665.  
  666. static void Transform_Superellipsoid(Object, Trans)
  667. OBJECT *Object;
  668. TRANSFORM *Trans;
  669. {
  670.   Compose_Transforms(((SUPERELLIPSOID *)Object)->Trans, Trans);
  671.  
  672.   Compute_Superellipsoid_BBox((SUPERELLIPSOID *)Object);
  673. }
  674.  
  675.  
  676.  
  677. /*****************************************************************************
  678. *
  679. * FUNCTION
  680. *
  681. *   Invert_Superellipsoid
  682. *
  683. * INPUT
  684. *
  685. *   Object - Object
  686. *   
  687. * OUTPUT
  688. *
  689. *   Object
  690. *   
  691. * RETURNS
  692. *   
  693. * AUTHOR
  694. *
  695. *   Dieter Bayer
  696. *   
  697. * DESCRIPTION
  698. *
  699. *   Invert a superellipsoid.
  700. *
  701. * CHANGES
  702. *
  703. *   Oct 1994 : Creation.
  704. *
  705. ******************************************************************************/
  706.  
  707. static void Invert_Superellipsoid(Object)
  708. OBJECT *Object;
  709. {
  710.   Invert_Flag(Object, INVERTED_FLAG);
  711. }
  712.  
  713.  
  714.  
  715. /*****************************************************************************
  716. *
  717. * FUNCTION
  718. *
  719. *   Create_Superellipsoid
  720. *
  721. * INPUT
  722. *   
  723. * OUTPUT
  724. *   
  725. * RETURNS
  726. *
  727. *   SUPERELLIPSOID * - new superellipsoid
  728. *   
  729. * AUTHOR
  730. *
  731. *   Dieter Bayer
  732. *   
  733. * DESCRIPTION
  734. *
  735. *   Create a new superellipsoid.
  736. *
  737. * CHANGES
  738. *
  739. *   Oct 1994 : Creation.
  740. *
  741. ******************************************************************************/
  742.  
  743. SUPERELLIPSOID *Create_Superellipsoid()
  744. {
  745.   SUPERELLIPSOID *New;
  746.  
  747.   New = (SUPERELLIPSOID *)POV_MALLOC(sizeof(SUPERELLIPSOID), "superellipsoid");
  748.  
  749.   INIT_OBJECT_FIELDS(New,SUPERELLIPSOID_OBJECT,&Superellipsoid_Methods)
  750.  
  751.   New->Trans = Create_Transform();
  752.  
  753.   Make_Vector(New->Power, 2.0, 2.0, 2.0);
  754.  
  755.   return(New);
  756. }
  757.  
  758.  
  759.  
  760. /*****************************************************************************
  761. *
  762. * FUNCTION
  763. *
  764. *   Copy_Superellipsoid
  765. *
  766. * INPUT
  767. *
  768. *   Object - Object
  769. *   
  770. * OUTPUT
  771. *   
  772. * RETURNS
  773. *
  774. *   void * - New superellipsoid
  775. *   
  776. * AUTHOR
  777. *
  778. *   Dieter Bayer
  779. *   
  780. * DESCRIPTION
  781. *
  782. *   Copy a superellipsoid structure.
  783. *
  784. * CHANGES
  785. *
  786. *   Oct 1994 : Creation.
  787. *
  788. ******************************************************************************/
  789.  
  790. static void *Copy_Superellipsoid(Object)
  791. OBJECT *Object;
  792. {
  793.   SUPERELLIPSOID *New, *Superellipsoid = (SUPERELLIPSOID *)Object;
  794.  
  795.   New = Create_Superellipsoid();
  796.  
  797.   /* Get rid of the transformation created in Create_Superellipsoid(). */
  798.  
  799.   Destroy_Transform(New->Trans);
  800.  
  801.   /* Copy superellipsoid. */
  802.  
  803.   *New = *Superellipsoid;
  804.  
  805.   New->Trans = Copy_Transform(Superellipsoid->Trans);
  806.  
  807.   return(New);
  808. }
  809.  
  810.  
  811.  
  812. /*****************************************************************************
  813. *
  814. * FUNCTION
  815. *
  816. *   Destroy_Superellipsoid
  817. *
  818. * INPUT
  819. *
  820. *   Object - Object
  821. *   
  822. * OUTPUT
  823. *
  824. *   Object
  825. *   
  826. * RETURNS
  827. *   
  828. * AUTHOR
  829. *
  830. *   Dieter Bayer
  831. *   
  832. * DESCRIPTION
  833. *
  834. *   Destroy a superellipsoid.
  835. *
  836. * CHANGES
  837. *
  838. *   Oct 1994 : Creation.
  839. *
  840. ******************************************************************************/
  841.  
  842. static void Destroy_Superellipsoid(Object)
  843. OBJECT *Object;
  844. {
  845.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  846.  
  847.   Destroy_Transform(Superellipsoid->Trans);
  848.  
  849.   POV_FREE (Object);
  850. }
  851.  
  852.  
  853.  
  854. /*****************************************************************************
  855. *
  856. * FUNCTION
  857. *
  858. *   Compute_Superellipsoid_BBox
  859. *
  860. * INPUT
  861. *
  862. *   Superellipsoid - Superellipsoid
  863. *   
  864. * OUTPUT
  865. *
  866. *   Superellipsoid
  867. *   
  868. * RETURNS
  869. *   
  870. * AUTHOR
  871. *
  872. *   Dieter Bayer
  873. *   
  874. * DESCRIPTION
  875. *
  876. *   Calculate the bounding box of a superellipsoid.
  877. *
  878. * CHANGES
  879. *
  880. *   Oct 1994 : Creation.
  881. *
  882. ******************************************************************************/
  883.  
  884. void Compute_Superellipsoid_BBox(Superellipsoid)
  885. SUPERELLIPSOID *Superellipsoid;
  886. {
  887.   Make_BBox(Superellipsoid->BBox, -1.0001, -1.0001, -1.0001, 2.0002, 2.0002, 2.0002);
  888.  
  889.   Recompute_BBox(&Superellipsoid->BBox, Superellipsoid->Trans);
  890. }
  891.  
  892.  
  893.  
  894. /*****************************************************************************
  895. *
  896. * FUNCTION
  897. *
  898. *   intersect_box
  899. *
  900. * INPUT
  901. *
  902. *   P, D       - Ray origin and direction
  903. *   dmin, dmax - Intersection depths
  904. *   
  905. * OUTPUT
  906. *
  907. *   dmin, dmax
  908. *   
  909. * RETURNS
  910. *
  911. *   int - TRUE, if hit
  912. *   
  913. * AUTHOR
  914. *
  915. *   Dieter Bayer
  916. *   
  917. * DESCRIPTION
  918. *
  919. *   Intersect a ray with an axis aligned unit box.
  920. *
  921. * CHANGES
  922. *
  923. *   Oct 1994 : Creation.
  924. *
  925. ******************************************************************************/
  926.  
  927. static int intersect_box(P, D, dmin, dmax)
  928. VECTOR P, D;
  929. DBL *dmin, *dmax;
  930. {
  931.   DBL tmin = 0.0, tmax = 0.0;
  932.  
  933.   /* Left/right. */
  934.  
  935.   if (fabs(D[X]) > EPSILON)
  936.   {
  937.     if (D[X] > EPSILON)
  938.     {
  939.       *dmin = (MIN_VALUE - P[X]) / D[X];
  940.  
  941.       *dmax = (MAX_VALUE - P[X]) / D[X];
  942.  
  943.       if (*dmax < EPSILON) return(FALSE);
  944.     }
  945.     else
  946.     {
  947.       *dmax = (MIN_VALUE - P[X]) / D[X];
  948.  
  949.       if (*dmax < EPSILON) return(FALSE);
  950.  
  951.       *dmin = (MAX_VALUE - P[X]) / D[X];
  952.     }
  953.  
  954.     if (*dmin > *dmax) return(FALSE);
  955.   }
  956.   else
  957.   {
  958.     if ((P[X] < MIN_VALUE) || (P[X] > MAX_VALUE))
  959.     {
  960.       return(FALSE);
  961.     }
  962.  
  963.     *dmin = -BOUND_HUGE;
  964.     *dmax =  BOUND_HUGE;
  965.   }
  966.  
  967.   /* Top/bottom. */
  968.  
  969.   if (fabs(D[Y]) > EPSILON)
  970.   {
  971.     if (D[Y] > EPSILON)
  972.     {
  973.       tmin = (MIN_VALUE - P[Y]) / D[Y];
  974.  
  975.       tmax = (MAX_VALUE - P[Y]) / D[Y];
  976.     }
  977.     else
  978.     {
  979.       tmax = (MIN_VALUE - P[Y]) / D[Y];
  980.  
  981.       tmin = (MAX_VALUE - P[Y]) / D[Y];
  982.     }
  983.  
  984.     if (tmax < *dmax)
  985.     {
  986.       if (tmax < EPSILON) return(FALSE);
  987.  
  988.       if (tmin > *dmin)
  989.       {
  990.         if (tmin > tmax) return(FALSE);
  991.  
  992.         *dmin = tmin;
  993.       }
  994.       else
  995.       {
  996.         if (*dmin > tmax) return(FALSE);
  997.       }
  998.  
  999.       *dmax = tmax;
  1000.     }
  1001.     else
  1002.     {
  1003.       if (tmin > *dmin)
  1004.       {
  1005.         if (tmin > *dmax) return(FALSE);
  1006.  
  1007.         *dmin = tmin;
  1008.       }
  1009.     }
  1010.   }
  1011.   else
  1012.   {
  1013.     if ((P[Y] < MIN_VALUE) || (P[Y] > MAX_VALUE))
  1014.     {
  1015.       return(FALSE);
  1016.     }
  1017.   }
  1018.  
  1019.   /* Front/back. */
  1020.  
  1021.   if (fabs(D[Z]) > EPSILON)
  1022.   {
  1023.     if (D[Z] > EPSILON)
  1024.     {
  1025.       tmin = (MIN_VALUE - P[Z]) / D[Z];
  1026.  
  1027.       tmax = (MAX_VALUE - P[Z]) / D[Z];
  1028.     }
  1029.     else
  1030.     {
  1031.       tmax = (MIN_VALUE - P[Z]) / D[Z];
  1032.  
  1033.       tmin = (MAX_VALUE - P[Z]) / D[Z];
  1034.     }
  1035.  
  1036.     if (tmax < *dmax)
  1037.     {
  1038.       if (tmax < EPSILON) return(FALSE);
  1039.  
  1040.       if (tmin > *dmin)
  1041.       {
  1042.         if (tmin > tmax) return(FALSE);
  1043.  
  1044.         *dmin = tmin;
  1045.       }
  1046.       else
  1047.       {
  1048.         if (*dmin > tmax) return(FALSE);
  1049.       }
  1050.  
  1051.       *dmax = tmax;
  1052.     }
  1053.     else
  1054.     {
  1055.       if (tmin > *dmin)
  1056.       {
  1057.         if (tmin > *dmax) return(FALSE);
  1058.  
  1059.         *dmin = tmin;
  1060.       }
  1061.     }
  1062.   }
  1063.   else
  1064.   {
  1065.     if ((P[Z] < MIN_VALUE) || (P[Z] > MAX_VALUE))
  1066.     {
  1067.       return(FALSE);
  1068.     }
  1069.   }
  1070.  
  1071.   *dmin = tmin;
  1072.   *dmax = tmax;
  1073.  
  1074.   return(TRUE);
  1075. }
  1076.  
  1077.  
  1078.  
  1079. /*****************************************************************************
  1080. *
  1081. * FUNCTION
  1082. *
  1083. *   evaluate_superellipsoid
  1084. *
  1085. * INPUT
  1086. *
  1087. *   P          - Point to evaluate
  1088. *   Superellipsoid - Pointer to superellipsoid
  1089. *   
  1090. * OUTPUT
  1091. *   
  1092. * RETURNS
  1093. *
  1094. *   DBL
  1095. *   
  1096. * AUTHOR
  1097. *
  1098. *   Dieter Bayer
  1099. *   
  1100. * DESCRIPTION
  1101. *
  1102. *   Get superellipsoid value in the given point.
  1103. *
  1104. * CHANGES
  1105. *
  1106. *   Oct 1994 : Creation.
  1107. *
  1108. ******************************************************************************/
  1109.  
  1110. static DBL evaluate_superellipsoid(P, Superellipsoid)
  1111. VECTOR P;
  1112. SUPERELLIPSOID *Superellipsoid;
  1113. {
  1114.   VECTOR V1;
  1115.  
  1116.   V1[X] = power(fabs(P[X]), Superellipsoid->Power[X]);
  1117.   V1[Y] = power(fabs(P[Y]), Superellipsoid->Power[X]);
  1118.   V1[Z] = power(fabs(P[Z]), Superellipsoid->Power[Z]);
  1119.  
  1120.   return(power(V1[X] + V1[Y], Superellipsoid->Power[Y]) + V1[Z] - 1.0);
  1121. }
  1122.  
  1123.  
  1124.  
  1125. /*****************************************************************************
  1126. *
  1127. * FUNCTION
  1128. *
  1129. *   power
  1130. *
  1131. * INPUT
  1132. *
  1133. *   x - Argument
  1134. *   e - Power
  1135. *   
  1136. * OUTPUT
  1137. *   
  1138. * RETURNS
  1139. *
  1140. *   DBL
  1141. *   
  1142. * AUTHOR
  1143. *
  1144. *   Dieter Bayer
  1145. *   
  1146. * DESCRIPTION
  1147. *
  1148. *   Raise x to the power of e.
  1149. *
  1150. * CHANGES
  1151. *
  1152. *   Oct 1994 : Creation.
  1153. *
  1154. ******************************************************************************/
  1155.  
  1156. static DBL power(x, e)
  1157. DBL x, e;
  1158. {
  1159.   register int i;
  1160.   register DBL b;
  1161.  
  1162.   b = x;
  1163.  
  1164.   i = (int)e;
  1165.  
  1166.   /* Test if we have an integer power. */
  1167.  
  1168.   if (e == (DBL)i)
  1169.   {
  1170.     switch (i)
  1171.     {
  1172.       case 0: return(1.0);
  1173.  
  1174.       case 1: return(b);
  1175.  
  1176.       case 2: return(Sqr(b));
  1177.  
  1178.       case 3: return(Sqr(b) * b);
  1179.  
  1180.       case 4: b *= b; return(Sqr(b));
  1181.  
  1182.       case 5: b *= b; return(Sqr(b) * x);
  1183.  
  1184.       case 6: b *= b; return(Sqr(b) * b);
  1185.  
  1186.       default: return(pow(x, e));
  1187.     }
  1188.   }
  1189.   else
  1190.   {
  1191.     return(pow(x, e));
  1192.   }
  1193. }
  1194.  
  1195.  
  1196.  
  1197. /*****************************************************************************
  1198. *
  1199. * FUNCTION
  1200. *
  1201. *   insert_hit
  1202. *
  1203. * INPUT
  1204. *
  1205. *   Object      - Object
  1206. *   Ray         - Ray
  1207. *   Depth       - Intersection depth
  1208. *   Depth_Stack - Intersection stack
  1209. *   
  1210. * OUTPUT
  1211. *
  1212. *   Depth_Stack
  1213. *   
  1214. * RETURNS
  1215. *
  1216. *   int - TRUE, if intersection is valid
  1217. *   
  1218. * AUTHOR
  1219. *
  1220. *   Dieter Bayer
  1221. *   
  1222. * DESCRIPTION
  1223. *
  1224. *   Push an intersection onto the depth stack if it is valid.
  1225. *
  1226. * CHANGES
  1227. *
  1228. *   Nov 1994 : Creation.
  1229. *
  1230. ******************************************************************************/
  1231.  
  1232. static int insert_hit(Object, Ray, Depth, Depth_Stack)
  1233. OBJECT *Object;
  1234. RAY *Ray;
  1235. DBL Depth;
  1236. ISTACK *Depth_Stack;
  1237. {
  1238.   VECTOR IPoint;
  1239.  
  1240.   if ((Depth > DEPTH_TOLERANCE) && (Depth < Max_Distance))
  1241.   {
  1242.     VEvaluateRay(IPoint, Ray->Initial, Depth, Ray->Direction);
  1243.  
  1244.     if (Point_In_Clip(IPoint, Object->Clip))
  1245.     {
  1246.       push_entry(Depth, IPoint, Object, Depth_Stack);
  1247.  
  1248.       return(TRUE);
  1249.     }
  1250.   }
  1251.  
  1252.   return(FALSE);
  1253. }
  1254.  
  1255.  
  1256.  
  1257. /*****************************************************************************
  1258. *
  1259. * FUNCTION
  1260. *
  1261. *   compdists
  1262. *
  1263. * INPUT
  1264. *   
  1265. * OUTPUT
  1266. *   
  1267. * RETURNS
  1268. *   
  1269. * AUTHOR
  1270. *
  1271. *   Alexander Enzmann
  1272. *   
  1273. * DESCRIPTION
  1274. *
  1275. *   Compare two slabs.
  1276. *
  1277. * CHANGES
  1278. *
  1279. *   Nov 1994 : Creation.
  1280. *
  1281. ******************************************************************************/
  1282.  
  1283. static int compdists(in_a, in_b)
  1284.  CONST void *in_a;
  1285.  CONST void *in_b;
  1286. {
  1287.   DBL a, b;
  1288.  
  1289.   a = *((DBL *)in_a);
  1290.   b = *((DBL *)in_b);
  1291.  
  1292.   if (a < b)
  1293.   {
  1294.     return(-1);
  1295.   }
  1296.  
  1297.   if (a == b)
  1298.   {
  1299.     return(0);
  1300.   }
  1301.   else
  1302.   {
  1303.     return(1);
  1304.   }
  1305. }
  1306.  
  1307.  
  1308.  
  1309. /*****************************************************************************
  1310. *
  1311. * FUNCTION
  1312. *
  1313. *   find_ray_plane_points
  1314. *
  1315. * INPUT
  1316. *   
  1317. * OUTPUT
  1318. *   
  1319. * RETURNS
  1320. *   
  1321. * AUTHOR
  1322. *
  1323. *   Alexander Enzmann
  1324. *   
  1325. * DESCRIPTION
  1326. *
  1327. *   Find all the places where the ray intersects the set of
  1328. *   subdividing planes through the superquadric.  Return the
  1329. *   number of valid hits (within the bounding box).
  1330. *
  1331. * CHANGES
  1332. *
  1333. *   Nov 1994 : Creation.
  1334. *
  1335. ******************************************************************************/
  1336.  
  1337. static int find_ray_plane_points(P, D, cnt, dists, mindist, maxdist)
  1338. VECTOR P, D;
  1339. int cnt;
  1340. DBL *dists, mindist, maxdist;
  1341. {
  1342.   int i;
  1343.   DBL t, d;
  1344.  
  1345.   /* Since min and max dist are the distance to two of the bounding planes
  1346.      we are considering, there is a high probablity of missing them due to
  1347.      round off error. Therefore we adjust min and max. */
  1348.  
  1349.   t = EPSILON * (maxdist - mindist);
  1350.  
  1351.   mindist -= t;
  1352.   maxdist += t;
  1353.  
  1354.   /* Check the sets of planes that cut apart the superquadric. */
  1355.  
  1356.   for (i = 0; i < PLANECOUNT; i++)
  1357.   {
  1358.     d = (D[0] * planes[i][0] + D[1] * planes[i][1] + D[2] * planes[i][2]);
  1359.  
  1360.     if (fabs(d) < EPSILON)
  1361.     {
  1362.       /* Can't possibly get a hit for this combination of ray and plane. */
  1363.  
  1364.       continue;
  1365.     }
  1366.  
  1367.     t = (planes[i][3] - (P[0] * planes[i][0] + P[1] * planes[i][1] + P[2] * planes[i][2])) / d;
  1368.  
  1369.     if ((t >= mindist) && (t <= maxdist))
  1370.     {
  1371.       dists[cnt++] = t;
  1372.     }
  1373.   }
  1374.  
  1375.   /* Sort the results for further processing. */
  1376.  
  1377.   QSORT((void *)(dists), (size_t)cnt, sizeof(DBL), compdists);
  1378.  
  1379.   return(cnt);
  1380. }
  1381.  
  1382.  
  1383.  
  1384. /*****************************************************************************
  1385. *
  1386. * FUNCTION
  1387. *
  1388. *   solve_hit1
  1389. *
  1390. * INPUT
  1391. *   
  1392. * OUTPUT
  1393. *   
  1394. * RETURNS
  1395. *   
  1396. * AUTHOR
  1397. *
  1398. *   Alexander Enzmann
  1399. *   
  1400. * DESCRIPTION
  1401. *
  1402. *   Home in on the root of a superquadric using a combination of
  1403. *   secant and bisection methods.  This routine requires that
  1404. *   the sign of the function be different at P0 and P1, it will
  1405. *   fail drastically if this isn't the case.
  1406. *
  1407. * CHANGES
  1408. *
  1409. *   Nov 1994 : Creation.
  1410. *
  1411. ******************************************************************************/
  1412.  
  1413. static void solve_hit1(Super, v0, tP0, v1, tP1, P)
  1414. SUPERELLIPSOID *Super;
  1415. DBL v0, v1;
  1416. VECTOR tP0, tP1, P;
  1417. {
  1418.   int i;
  1419.   DBL x, v2, v3;
  1420.   VECTOR P0, P1, P2, P3;
  1421.  
  1422.   Assign_Vector(P0, tP0);
  1423.   Assign_Vector(P1, tP1);
  1424.  
  1425.   /* The sign of v0 and v1 changes between P0 and P1, this
  1426.      means there is an intersection point in there somewhere. */
  1427.  
  1428.   for (i = 0; i < MAX_ITERATIONS; i++)
  1429.   {
  1430.     if (fabs(v0) < ZERO_TOLERANCE)
  1431.     {
  1432.       /* Near point is close enough to an intersection - just use it. */
  1433.  
  1434.       Assign_Vector(P, P0);
  1435.  
  1436.       break;
  1437.     }
  1438.  
  1439.     if (fabs(v1) < ZERO_TOLERANCE)
  1440.     {
  1441.       /* Far point is close enough to an intersection. */
  1442.  
  1443.       Assign_Vector(P, P1);
  1444.  
  1445.       break;
  1446.     }
  1447.  
  1448.     /* Look at the chord connecting P0 and P1. */
  1449.  
  1450.     /* Assume a line between the points. */
  1451.  
  1452.     x = fabs(v0) / fabs(v1 - v0);
  1453.  
  1454.     VSub(P2, P1, P0);
  1455.  
  1456.     VAddScaled(P2, P0, x, P2);
  1457.  
  1458.     v2 = evaluate_superellipsoid(P2, Super);
  1459.  
  1460.     /* Look at the midpoint between P0 and P1. */
  1461.  
  1462.     VSub(P3, P1, P0);
  1463.  
  1464.     VAddScaled(P3, P0, 0.5, P3);
  1465.  
  1466.     v3 = evaluate_superellipsoid(P3, Super);
  1467.  
  1468.     if (v0 * v2 > 0.0)
  1469.     {
  1470.       if (v1 * v2 > 0.0)
  1471.       {
  1472.         /* This should be impossible, since v0 and v1 were opposite signs,
  1473.            v2 must be either 0 or opposite in sign to either v0 or v1. */
  1474.  
  1475.         Error("internal failure in function solve_sq_hit1: %d, %g, %g, %g", i, v0, v1, v2);
  1476.       }
  1477.       else
  1478.       {
  1479.         if (v0 * v3 > 0.0)
  1480.         {
  1481.           if (x < 0.5)
  1482.           {
  1483.             Assign_Vector(P0, P3);
  1484.           }
  1485.           else
  1486.           {
  1487.             Assign_Vector(P0, P2);
  1488.           }
  1489.         }
  1490.         else
  1491.         {
  1492.           /* We can move both ends. */
  1493.  
  1494.           Assign_Vector(P0, P2);
  1495.  
  1496.           Assign_Vector(P1, P3);
  1497.         }
  1498.       }
  1499.     }
  1500.     else
  1501.     {
  1502.       if (v0 * v3 > 0.0)
  1503.       {
  1504.         /* We can move both ends. */
  1505.  
  1506.         Assign_Vector(P0, P3);
  1507.  
  1508.         Assign_Vector(P1, P2);
  1509.       }
  1510.       else
  1511.       {
  1512.         if (x < 0.5)
  1513.         {
  1514.           Assign_Vector(P1, P2);
  1515.         }
  1516.         else
  1517.         {
  1518.           Assign_Vector(P1, P3);
  1519.         }
  1520.       }
  1521.     }
  1522.   }
  1523.  
  1524.   if (i == MAX_ITERATIONS)
  1525.   {
  1526.     /* The loop never quite closed in on the result - just use the point
  1527.        closest to zero.  This really shouldn't happen since the max number
  1528.        of iterations is enough to converge with straight bisection.  */
  1529.  
  1530.     if (fabs(v0) < fabs(v1))
  1531.     {
  1532.       Assign_Vector(P, P0);
  1533.     }
  1534.     else
  1535.     {
  1536.       Assign_Vector(P, P1);
  1537.     }
  1538.   }
  1539. }
  1540.  
  1541.  
  1542.  
  1543.  
  1544. /*****************************************************************************
  1545. *
  1546. * FUNCTION
  1547. *
  1548. *   check_hit2
  1549. *
  1550. * INPUT
  1551. *   
  1552. * OUTPUT
  1553. *   
  1554. * RETURNS
  1555. *   
  1556. * AUTHOR
  1557. *
  1558. *   Alexander Enzmann
  1559. *   
  1560. * DESCRIPTION
  1561. *
  1562. *   Try to find the root of a superquadric using Newtons method.
  1563. *
  1564. * CHANGES
  1565. *
  1566. *   Nov 1994 : Creation.
  1567. *
  1568. ******************************************************************************/
  1569.  
  1570. static int check_hit2(Super, P, D, t0, P0, v0, t1, t, Q)
  1571. SUPERELLIPSOID *Super;
  1572. DBL t0, t1, v0, *t;
  1573. VECTOR P, D, P0, Q;
  1574. {
  1575.   int i;
  1576.   DBL dt0, dt1, v1, deltat, maxdelta;
  1577.   VECTOR P1;
  1578.  
  1579.   dt0 = t0;
  1580.   dt1 = t0 + 0.0001 * (t1 - t0);
  1581.  
  1582.   maxdelta = t1 - t0;
  1583.  
  1584.   for (i = 0; (dt0 < t1) && (i < MAX_ITERATIONS); i++)
  1585.   {
  1586.     VEvaluateRay(P1, P, dt1, D)
  1587.  
  1588.     v1 = evaluate_superellipsoid(P1, Super);
  1589.  
  1590.     if (v0 * v1 < 0)
  1591.     {
  1592.       /* Found a crossing point, go back and use normal root solving. */
  1593.  
  1594.       solve_hit1(Super, v0, P0, v1, P1, Q);
  1595.  
  1596.       VSub(P0, Q, P);
  1597.  
  1598.       VLength(*t, P0);
  1599.  
  1600.       return(TRUE);
  1601.     }
  1602.     else
  1603.     {
  1604.       if (fabs(v1) < ZERO_TOLERANCE)
  1605.       {
  1606.          VEvaluateRay(Q, P, dt1, D)
  1607.  
  1608.          *t = dt1;
  1609.  
  1610.          return(TRUE);
  1611.       }
  1612.       else
  1613.       {
  1614.         if (((v0 > 0.0) && (v1 > v0)) || ((v0 < 0.0) && (v1 < v0)))
  1615.         {
  1616.           /* We definitely failed */
  1617.  
  1618.           break;
  1619.         }
  1620.         else
  1621.         {
  1622.           if (v1 == v0)
  1623.           {
  1624.             break;
  1625.           }
  1626.           else
  1627.           {
  1628.             deltat = v1 * (dt1 - dt0) / (v1 - v0);
  1629.           }
  1630.         }
  1631.       }
  1632.     }
  1633.  
  1634.     if (fabs(deltat) > maxdelta)
  1635.     {
  1636.        break;
  1637.     }
  1638.  
  1639.     v0 = v1;
  1640.  
  1641.     dt0 = dt1;
  1642.  
  1643.     dt1 -= deltat;
  1644.   }
  1645.  
  1646.   return(FALSE);
  1647. }
  1648.